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Abstract The Discrete Particle Method (DPM) is used to 
model granular flows down an inclined chute. We observe 
three major regimes: static piles, steady uniform flows and 
accelerating flows. For flows over a smooth base, other (quasi- 
steady) regimes are observed where the flow is either highly 
energetic and strongly layered in depth for small inclina- 
tions, or non-uniform and oscillating for larger inclinations. 

For steady uniform flows, depth profiles of density, ve- 
locity and stress have been obtained using an improved coarse- 
graining method, which allows accurate statistics even at the 
base of the flow. A shallow-layer model for granular flows is 
completed with macro-scale closure relations obtained from 
micro-scale DPM simulations of steady flows. We thus ob- 
tain relations for the effective basal friction, shape factor, 
mean density, and the normal stress anisotropy as functions 
of layer thickness, flow velocity and basal roughness. For 
collisional flows, the functional dependencies are well de- 
termined and have been obtained. 

Keywords Discrete Particle Method • Coarse graining • 
Granular chute flow • Depth-averaging • Shallow-layer 
equations 



1 General introduction 

1.1 Background 

Granular avalanche flows are common to natural environ- 
ments and industry. They occur across many orders of mag- 
nitude. Examples range from rock slides, containing upwards 
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of lOOOm^ of material; to the flow of sinter, pellets and coke 
into a blast furnace for iron-ore melting; down to the flow of 
sand in an hour-glass. The dynamics of these flows are in- 
fluenced by many factors such as: polydisperity; variations 
in density; non-uniform shape; complex basal topography; 
surface contact properties; coexistence of static, steady and 
accelerating material; and, flow obstacles and constrictions. 

Discrete Particle Methods (DPMs) are an extremely pow- 
erful way to investigate the effects of these and other factors. 
With the rapid recent improvement in computational power 
the full simulation of the flow in a small hour glass of mil- 
lions of particles is now feasible. However, complete DPM 
simulations of large-scale geophysical mass flow will, prob- 
ably, never be possible. 

One primary goal of the present research is to simu- 
late large scale and complex industrial flows using granular 
shallow-layer equations. In this paper we will take the first 
step of using the DPM ^ CS79llSEG+01IISGPL02IISLG03l 
iLudOSI to simulate small granular flows of mono-dispersed 
spherical particles in steady flow situations. We will use a 
refined and novel analysis to investigate three particular as- 
pects of shallow chute flows: i) how to obtain meaningful 
macro-scale fields from the DPM simulation, ii) how to asses 
the flow dependence on the basal roughness, and Hi) how to 
validate the assumptions made in depth-averaged theory. 

The DPM simulations presented here will enable the con- 
struction of the mapping between the micro-scale and macro- 
scale variables and functions, thus enabling construction of a 
closed set of continuum equations. These mappings (closure 
relations) can then be used in continuum shallow-layer mod- 
els and compared with full DPM simulations (DPMs). For 
certain situations, precomputed closures should work; but, 
in very complicated situations the pre-established relations 
may fail. Heterogeneous, multi-scale modelling (HMM) is 
then an alternative liWEL+07j . in which the local consiti- 
tute relations are directly used in the continuum model. In 
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HMM, continuum and micro-scale models are dynamically 
coupled with a two-way communication between the differ- 
ent models in selective regions in both space and time, thus 
reducing computational expense and allowing simulation of 
complex granular flows. 

1 .2 Shallow-layer models 

Shallow granular continuum models are often used to sim- 
ulate geophysical mass flows, including snow avalanches 
IICGJ07I . dense pyroclastic flows, debris flows OPTO II . block 
and ash flows flPPP+OSl and lahars IIWSS0 8 I. Such shaflow- 
layer models involve approximations reducing the proper- 
ties of a huge number of individual particles to a handful 
of averaged quantities. Originally these models were de- 
rived from the general continuum incompressible mass and 
momentum equations, using the long-wave approximation 
llSH89llHSSN93IIGWH99llWGH99IIGTN03llBTT2l for shal- 
low variations in the flow height and basal topography. De- 
spite the massive reduction in degrees of freedom made, 
shallow-layer models tend to be surprisingly accurate, and 
are thus an effective tool in modelling geophysical flows. 
Moreover, they are now used as a geological risk assessment 
and hazard planning tool llDPP+081 . In addition to these 
geological applications, model applications involve small- 
scale laboratory chute flows containing obstacles IIGTN03I . 
wedges IHH05l lGC07l and conti-actions IVATK+071. show- 
ing good quantitative agreement between theory and exper- 
iment. 

In fluid dynamics, the Navier-Stokes equations are es- 
tablished with full constitutive equations. Nonetheless, the 
shallow-layer equations or Saint- Venant equations are often 
used in large scale situations where it is unpractical to solve 
the full Navier-Stokes equations. Our present aim is to di- 
rectly investigate the validity of the assumptions of granu- 
lar shallow-layer models first from discrete particle simula- 
tions, before obtaining fully 3D 'kinetic theory' -style consti- 
tutive relations and simplifying these via the depth-integration 
process. A discussion of the full three-dimensional prop- 
erties of our particle simulations will be undertaken later 
Here we restrict attention to the closures required for two- 
dimensional shallow granular flow equations. 

A key difference between shallow-layer fluid models and 
granular ones is the appearance of a basal friction coeffi- 
cient, jU, being the ratio of the shear to normal traction at the 
base. In early models, a dry Coulomb-like friction law was 
used IISH89I . It implies ji to be constant, given by the tan- 
gent of the friction angle between the material and the base, 
5, i.e., ji — tan5. As a consequence constant uniform flow 
is only possible in such a model at that angle 5, independent 
of height. There is a considerable amount of experimental 
evidence, e.g., P DD99IIGDR04J . that suggests that such a 
simple Coulomb law does not hold on rough beds or for 



moderate inclination angles. Furthermore, detailed experi- 
mental investigations using glass beads IPou99ll lead to an 
improved empirical 'Pouliquen' friction law characterised 
by two angles: the angle at which the material comes to rest, 
5i, where friction dominates over gravity and the angle, ^2, 
above which gravity dominates over friction and the material 
accelerates. Between these two angles steady flow is possi- 
ble, and in the limit 5i — > ^2 = 5 the original Coulomb style 
model is recovered. 

Since its formulation a lot of work has been performed 
on extending and understanding this Pouliquen law. The orig- 
inal law was obtained by retarding flowing material and mea- 
suring the angle at which the material stopped as a function 
of height hstup (6), or equivalently, by inverting this relation, 
Ostop{h)- For most materials, granular included, a greater 
angle is required to initiate stationary than to retard flow- 
ing material. Pouliquen and Forterre OPF02I . by measuring 
the angle required to start motion, measured Ostart (h), i-e., 
the friction law for initially stationary material. As expected 
Ostart was greater than Ostop and this information was used to 
extend the friction law to all values of the height and veloc- 
ity within the steady regime. Borzsonyi & Ecke performed 
a series of additional experiments: firstly, in |BE06| they 
looked at higher angles were the mean flow rates are close to 
the terminal velocity of a single particle, and found regions 
were the Pouliquen law is not valid and the Froude num- 
ber becomes inversely proportional to the height, as opposed 
to the linear relationship observed for steady flow. Borz- 
sonyi and Ecke, and Pouliquen and Forterre IIBE07IIFP03I 
have all worked on extending the original law to be valid 
for more complicated non-spherical materials like sand and 
metallic materials. Also, the effect of basal surface rough- 
ness has been systemically studied in IIGTDD031 by vary- 
ing the size of both the free flow and fixed basal particles. 
For convenience, we define X to be the size ratio of the 
fixed and the free particles. They observed a peak in rough- 
ness at a certain diameter ratio, A( , which depends on the 
compactness of the basal layer. Measured values of Xc in 
I GTDDQ3 I ranged between 1 and 3 for a monolayer of fixed 
particles. For fixed particles with smaller size and X < Xc, 
the range of angles where steady flow was observed de- 
creased, and eventually the steady flow regime completely 
vanished, i.e., ^2 — 5i— >OasA— >0 (yielding Coulomb type 
behaviour). For smaller flow particle diameters, i.e., with 
A > Xc, there was also a reduction in friction, but weaker 
than in the small X case. For much larger X, the friction sat- 
urated to a constant value, which they contributed to free 
particles that filled the holes in the basal surface and effec- 
tively created a stable basal surface of free particles. In a 
later publication MGDDT07I . they extended this investiga- 
tion to flows containing two particle sizes. 

Louge and Keast ILKOll modified the kinetic theory pre- 
sented in |,Jen93J by modelling enduring contacts via a fric- 
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tional rate-independent stress component in order to obtain 
steady flow on flat frictional inclines. This work was later 
extended to bumpy inclines IILou03l . Jenkins [Jen06l took 
a different approach and theoretically formulated a pheno- 
menological modification of granular kinetic theory to ac- 
count for enduring particle contacts. His idea is that endur- 
ing contacts between grains, forced by the shearing, reduce 
the collision rate of dissipation. Therefore a modification 
to the dissipation is introduced, which does not affect the 
stress. It leads to a law very similar to the one experimentally 
obtained by Pouliquen. He extended the theory in IIJen07l 
to very dissipative frictional particles, with a coefficient of 
restitution less than 0.7. Later, a detailed comparison with 
new experiments was performed, showing agreement with 
flows of low inclination [ JBIO]. 

Silbert et al [SEG+01| used DPMs to simulated chute 
flow of cohesionless particles. They found that a steady-state 
flow regime exists over a wide range of inclination angles, 
heights and interaction parameters, in confirmation of the 
experiments of Pouliquen IIPou99l . They found for steady- 
state flows that the volume fraction is constant throughout 
the flow, in agreement with the assumptions of shallow-layer 
theory IISH89I . They also observed that the shear stress is 
proportional to the square of the shear and the flow veloc- 
ity scales with the height to the power 3/2. This result co- 
incides with Bagnold's analysis of dilute binary collisions 
flows | Bag54[ . They also observed small systematic devia- 
tions from isotropic stress, which shows a deviation from 
fluid-like behaviour. However, normal stresses do not ap- 
proach a Coulomb-yield criterion structure at the angle of re- 
pose except near the surface, hinting that the failure of flow 
starts near the surface. In I SGPL02I . they looked at the effect 
of different basal types and found that for an ordered chute 
base the steady state regime splits into three distinct flow 
regimes: at smaller angles, the flowing system self-organizes 
into a state of low-dissipation flow consisting of in-plane 
ordering in the bulk; at higher angles, a high-dissipation 
regime similar to that for a rough base but with considerable 
slip at the bottom is observed; and, between these two sub- 
regions they observe a transitional flow regime characterized 
by large oscillations in the bulk averaged kinetic energy due 
to the spontaneous ordering and disordering of the system as 
a function of time. Finally, I1SLG03II investigated the initia- 
tion and cessation of granular chute flow more carefully and 
computed both Qstop and Qaun ■ For inclinations ^ Qstop 
they observed a Bagnold rheology, for ^ Q^top a linear pro- 
file, and for Q « Qstop avalanching behaviour 



1.3 Overview of this study 

Our present research is novel on the following three counts: 
Firstly, we compute more meaningful macro-scale fields 
from the DPM simulations than before by carefully choos- 



ing the coarse graining function. In order to homogenize the 
DPM data, the micro-scale fields need to be coarse-grained 
to obtain macroscopic fields. Coarse-grained micro-scale fields 
of density, momentum and stress were derived directly from 
the mass and momentum balance equations by Goldhirsch 
MGollOl . The quality of the statistics involved depends on the 
coarse graining width w, which defines the amount of spatial 
smoothing. For small coarse-graining width w, micro-scale 
variations remain visible, while for large w these smooth out 
in the macro-scale gradients. Since one of the objectives is to 
obtain the value of at the base, we use a novel adaptation 
of Goldhirsch' statistics |GollO| near boundaries. This is a 
non-trivial issue neglected in the literature; often continuum 
fields are simply not computed within a distance w of the 
boundary. Our new approach [ WTLB 1 1 II is consistent with 
the continuum equations everywhere, enabling the construc- 
tion of continuum fields within one course-graining width of 
the boundary. 

Secondly, we follow the approach of IIGTDD03I and vary 
the basal particle diameter to achieve different basal condi- 
tions. For particles with smaller basal than flowing diam- 
eter, A < 1, the flow becomes more energetic and oscilla- 
tory behaviour is observed. This phenomena has previously 
been reported in IISGPL02L but was achieved by changing 
the basal particles to a more regular, grid-like configuration. 
By investigating flow over fixed particles of different size 
than the free ones, we are able to quantify the roughness and 
numerically investigate the transition from rough to smooth 
surfaces. For smoother surfaces, we show that the parameter 
space can be split into to two types of steady flow, and we 
obtain a general friction law. 

Finally, we test the assumptions made in depth-averaged 
theory and determine the required closure laws. For shallow 
granular flows, the flow can be described by depth-averaged 
mass and momentum-balance equations [GTN031. Solving 
the depth-averaged equations requires a constitutive relation 
for the basal friction, a way to account for mean density 
variations, the shape of the velocity profile and the pressure 
anisotropy. We extract such data from DPMs obtained for 
steady uniform flows, and establish a novel, extended set of 
closure equations. Also, the depth-averaged equations are 
obtained under the assumptions that a) the density is con- 
stant in space and time and does not vary through the flow; 
b) the ratio between mean squared velocity and the squared 
mean velocity is known; c) the downward normal stress is 
lithostatic, i.e., balances the gravitational forces acting on 
the flow; and, d) the ratio between the normal stresses is 
known. Gray et al. IIGTN03II assumed the latter ratio to be 
one. The depth profiles of these quantities are discussed by 
Silbert et al. in ||SEG+01|ISGPL02IISLG03II for steady flow. 
We extend these measurements to validate our DPM simula- 
tions, using our superior statistical procedure. Hence, we es- 
tablish the range in which the shallow-layer model is valid. 
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and the closure relations required for shallow-layer contin- 
uum simulations, for a much wider range of flow regimes 
than had been considered before. 



1.4 OutUne 

We introduce the force model used in the DPM in Sect. |2l 
and the statistical method used to obtain macroscopic den- 
sity, velocity and stress profiles in Sect. [3] In Sect.|4l we dis- 
cuss the continuum shallow-layer equations for modelling 
granular flow including some macro-scale closures. The set 
up of the simulations is discussed in Sect.|5] and the steady- 
state regime is mapped for flows over a rough basal sur- 
face in Sect. |6] Depth profiles of the flow are introduced in 
Sect. |2l which are then used to characterize the steady flow 
over smoother surfaces in Sect. [8] Finally, the closure rela- 
tions for the shallow-layer model are established in Sect. |9] 
before we conclude in Sect. [10] 



2 Contact law description 

A Discrete Particle Method (DPM) is used to perform the 
simulation of a collection of identical granular particles. Bound- 
aries are created by special fixed particles, which generally 
will have different properties than the flow particles. Parti- 
cles interact by the standard spring-dashpot interaction model 
IICS79IIWB86llLud 081, in which it is assumed that particles 
are spherical and soft, and that pairs have a single contact 
point. 

Each particle / has a diameter dj, density p,-, position r,, 
velocity v,- and angular velocity a;, . For pairs of two particles 
we define the relative distance vector rij = r, — r,-, 
their separation rtj = |r,, |, the unit normal n,,- = rjj/rij and 
the relative velocity v,y = v; — v/. Two particles are in con- 
tact if their overlap. 



5^j = max{0,{di + dj)/2-rij), 



(1) 



is positive. A single contact point c at the centre of the over- 
lap is assumed, which is a valid assumption as long as the 
overlap is small. For our simulations the overlap between 
two particles is always below 1% of the particle radius; the 
simple contact model we employ thus captures the key pro- 
cess as there are no multiple contact points. 

The force acting on particle / is a combination of the 
pairwise interaction of two particles. The force f,j represents 
the force on particle ; from the interaction with particle j and 
can be decomposed into a normal and a tangential compo- 
nent. 



f f "■ + 1' ■ 

•i; — + 



(2) 



We assume that the particles are viscoelastic; therefore, 
they experience elastic as well as dissipative forces in both 



normal and tangential directions. The normal force is mod- 
elled as a spring-dashpot with a linear elastic and a linear 
dissipative contribution. 



f" — k"8 "n - -v"\"- 



(3) 



with spring constant k", damping coefficient y" and the nor- 
mal relative velocity component. 



(4) 



For a central collision, no tangential forces are present, and 
the collision time tc between two particles can be calculated 

as 



tc = %l\ 



k" ( 7" 



mij \2mij 



(5) 



with the reduced mass m,/ — mimj/{mi + mj). The normal 
restitution coefficient (ratio of relative normal speed after 
and before collision) is calculated as 



r,. = exp(-f,.77(2m,;))- 



(6) 



We also assume a linear elastic and a linear dissipative 
force in the tangential direction. 



f 

ij 



7% 



(7) 



with spring constant k' , damping coefficient 7', elastic tan- 
gential displacement 5'^ (which is explained later), and the 
relative tangential velocity at the contact point. 



(8) 



with b,j — — [idi — 5,")/2)n,y the branch vector from point 
/ to the contact point; for equal size particles b,; = — r,j/2. 

The elastic tangential force is modelled and derives from 
the roughness of the particle surface. Near the contact point, 
small bumps on a real particle would stick to each other, due 
to the normal force pressing them together, and elongate in 
the tangential direction resulting in an elastic force propor- 
tional to the elastic tangential displacement. It is defined to 
be zero at the initial time of contact, and its rate of change is 
given by 



d8'- id'--- 

IJ _ V? _ ^ 'J 

ij 



dt 



(9) 



where the first term is the relative tangential velocity at the 
contact point, and the second term ensures that SJy remains 
normal to n,;. The second term is always orthogonal to the 
spring direction and, hence, does not affect the rate of change 
of the spring length: it simply rotates it, thus keeping it tan- 
gential. 

When the tangential to normal force ratio becomes larger 
than the particle contact friction coefficient, /i^ , for a real 
particle the bumps would slip against each other and their 
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elongation is shortened until the bumps can stick to each 
other again. This is modelled by a static yield criterion, trun- 
cating the magnitude of 5\j as necessary to satisfy |f^| < 
IJ-c\f'ij\ - Thus, the contact surfaces are treated as stuck while 
\flj\ < IJ'c\f"j\ and as slipping otherwise, when the yield cri- 
terion is satisfiecQ. 

The total force on particle / is a combination of contact 
forces f,, with other particles and external forces such as 
gravity g. The resulting force f, and torque q, acting on par- 
ticle i Sire 

N N 

f,=g+ ™d L bii^f'J- (10) 

./=1../V'' j=lJ^i 

Finally, using these expressions we arrive at Newton's equa- 
tions of motion for the translational and rotational degrees 
of freedom, 

d^r, d 
Mi—^^fi, and /,— a;,= q;, (11) 
at^ at 

with m, the mass and /,■ the inertia of particle /. We integrate 
(fTTI) forward using Velocity- Verlet [AT93|, formally second 
order in time, with an adequate time step of At = /50. The 
collision time is given by (|5]l, while (|9]l is integrated using 
first-order forward Euler 

Hereafter, we distinguish between identical free flowing 
and identical fixed basal particles. Base particles are mod- 
elled as having an infinite mass and are unaffected by body 
forces: they do not move. This leaves two distinct types of 
collision: flow-flow, and flow-base. Model parameters for 
each of these collision types are changed independently. 

3 Statistics 

3.1 Coarse-graining 

The main aims of this paper are to use discrete particle sim- 
ulations to both confirm the assumptions of and provide clo- 
sure rules required for the depth-averaged shallow-water equa- 
tions. Hence, continuum fields have to be extracted from 
the discrete particle data. There are many papers in the lit- 
erature on how to go from the discrete to the continuum, 
binning micro-scale fields into small volumes [IK50 SH82j 
ILud04yLud09llLLV+0li . averaging along planes LTEDga, 
or coarse graining spatially and temporally OBab97IISA04l 
IGollOL Here, this will be achieved using a coarse-graining 
approach first described in |Bab97| based on a later deriva- 
tion of Goldhirsch IGollOI . extended further by us IIWTLBllI 

' Meant for review stage only. It should be noted that in the absence 
of dissipative forces and slipping, the system can be described as an 
Hamiltonian system: see Appendix|A] Appendix|B]contains details on 
the tangential displacement. A pseudocode of the tangential force cal- 
culation is provided in Appendix ICl 



to account for boundary forces due to the presence of the 
base. 

The method has several advantages over other methods 
because: (i) the fields produced automaticaUy satisfy the equa- 
tions of continuum mechanics, even near the flow base; (ii) 
it is neither assumed that the particles are rigid nor spheri- 
cal; and, (iii) the results are even valid for single particles as 
no averaging over groups of particles is required. The only 
assumptions are that each particle pair has a single point of 
contact {i.e., the particle shapes are convex), the contact area 
can be replaced by a contact point {i.e., the particles are not 
too soft), and that collisions are not instantaneous. 

3.2 Mass and momentum balance 
3.2.1 Notation and basic ideas 

Vectorial and tensorial components are denoted by Greek 
letters in order to distinguish them from the Latin particle 
indices i,j. Bold vector notation will be used when conve- 
nient. 

Assume a system given by Np flowing particles and A^^ 
fixed particles. Since we are interested in the flow, we will 
calculate macroscopic fields pertaining to the flowing parti- 
cles only. From statistical mechanics, the microscopic mass 
density of the flow, p™^'', at a point r at time t is defined by 

p'^'%r,t) = Y^mi5{r-ri{t)), (12) 

1=1 

where 5{r) is the Dirac delta function and m, is the mass 
of particle /. The following definition of the macroscopic 
density of the flow is used 

p{r,t) = Y,miW{r^ri{t)), (13) 

!=1 

thus replacing the Dirac delta function in (fT2l) by an in- 
tegrable 'coarse-graining' function W whose integral over 
space is unity. We will take the coarse-graining function to 
be a Gaussian 

with width or variance w. Other choices of the coarse-graining 
function are possible, but the Gaussian has the advantage 
that it produces smooth fields and the required integrals can 
be analyzed exactly. According to [GollOl, the answer de- 
pends only weakly on the choice of function, and the width 
w is the key parameter. 

It is clear that as w ^ the macroscopic density defined 
in ( fT4b reduces to the one in ( fT3l l. The coarse-graining func- 
tion can also be seen as a convolution integral between the 
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micro and macro definitions, 

p{r,t) ^ J W{r -r')p'^^''{r' ,t)dr' . 

3.2.2 Mass balance 



(15) 



Next we will consider how to obtain the other fields of in- 
terest: the momentum vector field and the stress tensor. As 
stated in Sect. 13. li the macroscopic variables will be defined 
in a way compatible with the continuum conservation laws. 
The coarse grained momentum density p{r,t) is defined 

by 



i=l 



(16) 



where the vta's are the velocity components of particle 
The macroscopic velocity field V{r,t) is then defined as 
the ratio of momentum and density fields, 



Vair,t):^pair,t)/p{r,t). 



(17) 



It is straightforward to confirm that equations (T3[ and ( fTSI l 
lead to the continuity equation 



dp , dpa 



0, 



(18) 



dt dra 

with the Einstein summation convention for Greek letters. 
3.2.3 Momentum balance 

Finally, we will consider the momentum conservation equa- 
tion with the aim of establishing the macroscopic stress field. 
In general, the desired momentum balance equations are writ- 
ten as. 



dpa 
dt 



da, 



'^[pVaVp]+^^+P8a 
dra dra 



(19) 



where is the stress tensor, and ga a component of the 
acceleration vector of gravity. 

Expressions ( fTSI l and ([TtI i for the momentum p and the 
velocity V have already been defined. The next step is to 
compute their temporal and spatial derivatives, respectively, 
and reach closure. Taking the time derivative of (fTST i gives 

^ ^^imv.a{tW{r~r,{t)) 



Np Np ^ 

^ mmaWir - r.) + miVia-^W{r - n). (20) 

/= 1 /= 1 



Using (fTTT i. the first term in (l20l l can be expressed as 

Np Np 

A„ = £ m,v,„^(r - n) = Y^fia^ir - r,). (21) 

i=l i=l 



In the simulations presented later the force on each par- 
ticle contains three contributions: particle-particle interac- 
tions, particle-base interactions, and the gravitational body 
force. Hence, 



Np Nh 

fia= Y fUa + Yfika+'^ygcc^ 

j=\,j^i k=\ 



(22) 



where fij is the interaction force between particle / and j, 
and f°^. the interaction between particle / and base particle k, 
or base wall if the base is flat. Therefore, we rework (1211 1 as 

Np Np Np N, 

i=ij=l,Mi !=U=1 i=l 

where Wi — W{r ~ ri). The last term in (l23T l can be simpli- 
fied to pga by using ( fT3l l. From Newton's third law, the con- 
tact forces are equal and opposite, such that fij = — /j,.Hence, 



i=ij=ij^i i=ij=i,i^j i^ij^ii^j 

(24) 

where in the first step we interchanged the dummy summa- 
tion indices. It follows from (l24b that ( l23T l can be written 

as 



I Np Np 



Np Nt 



ija Vi - rFj) 



i=\k=\ 



Np Nt 



-EE fija m - #:,) + E E /;'«^ + psa- (25) 

1=1, /■=!■+ 1 i=lk=\ 

Next, we will write Aa as the divergence of a tensor in 
order to obtain a formula for the stress tensor. The following 
identity holds for any smooth function W 

"1 d 



I —W{r-ri + srij)ds 
Jo ds 



^ drp Jo 



r-ri + srij)ds, 



(26) 



where Vij = r, — rj; we used the chain rule and differentia- 
tion to the full argument of per component. 

The next step extends Goldhirsch' analysis near a bound- 
ary. To obtain a similar expression for the interaction with 
base particles, we write 



Wi^ —^{r-ri + srik)ds 
Jo ds 



d 
dr^ 



/ W{r-ri + srik)ds, 
Jo 



(27) 
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which holds because #; decays towards infinity. Substitut- 
ing identities (|26ll, ^} and into (|25l) leads to 

(28) 



From BGollOI . it follows that the second term in ( l20l i can be 
expressed as follows 



E '"'^'a ^ #' (r - r,) = - 



pVaV/3+E"^'^'^«^';/3^ 



(29) 



where Vf is the fluctuating velocity of particle /, with com- 
ponents given by 

y;„(r,r) = y,„(f)-y„(r,0. (30) 
Substituting (l28T l and (|29] | into momentum balance (T9[ yields 



da, 



a/3 



E E /'■;«'"'■, 

i=l j=i+l 



UP 



{r — ri + srij)ds 



i=U=l 



(31) 



Therefore the stress is given by 

^«i5 = ~E E fijanip / 5r(r-ri + ir,y)di 

1=1 ;=!+l ■'0 



i=U=l 



(32) 



The kinetic component of the stress tensor as well as the 
contact stress coming from normal forces are symmetric; 
only the contact stress from tangential forces can contribute 
to the antisymmetric part of the stress tensor. In our simula- 
tions the tangential forces contribute less than 5% to the total 
stress in the system, such that the stress is almost symmetric. 

Equation ( |32] | differs from the results of IGollO l by an 
additional term that accounts for the stress created by the 
presence of the base. The contribution to the stress from the 
interaction of two flow particles /, /' is spatially distributed 
along the contact line from r; to r^, while the contribution 
from the interaction of particles ; with a fixed particle k is 
distributed along the line from r,- to r<-, extending further 
beyond r^.. We explain the situation as follows, see Fig. [1] 
Stress and density profiles are calculated using ( fTSl l and ( |32] | 




-2 



P 
a 



P 
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Fig. 1: Stress and density profiles are shown for two one- 
dimensional two-particle systems, each with two particles 
of unit mass at positions jc = ±1, and repelling each other 
(so with d>2 for our granular case). In the top figure, both 
particles are flowing, while in the bottom figure the left par- 
ticle is fixed and the right one flowing. 



for two one-dimensional two-particle systems, each with two 
particles of unit mass at positions jc = ±1, repelling each 
other with a force |/| = 1 and with w = 0.2. In the top figure, 
both particles belong to the flowing species, so the density is 
distributed around the particles' center of mass and the stress 
along the contact line. In the bottom figure, the left particle 
is a fixed base particle and the right particle is a free flowing 
one, so density is distributed around the flowing particle's 
center of mass and the stress along the line extending from 
the flowing particle to negative infinity. 

The strength of this statistical method is that the spa- 
tial coarse graining satisfies the mass and momentum bal- 
ance equations exactly at any given time, irrespective of the 
choice of the coarse graining function. Further details about 
the accuracy of the stress definition (|32] | are discussed in 
BWTLBl II . The expression for the energy is also not treated 
in this pubUcation, see MBab97l . 



4 Mathematical background 

In this section, we briefly outhne the existing knowledge to 
develop a continuum shallow-layer theory for granular flow. 



4.1 Shallow-layer model 

Shallow-layer models have been shown to be an effective 
tool in modelling many geophysical mass flows. Early ava- 
lanche models were formulated by adding gravitational ac- 
celeration and Coulomb basal friction to shallow-layer mod- 
els IIGEY67IIKE73I . Similar dry granular models have been 
derived using the long-wave approximation IISH89IIHSSN93I 
|Ive97llGWH99 WGH99 J for shaUow variations in flie flow 
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height and slope topography and included a Mohr-Coulomb 
rheology by the use of an earth pressure coefficient. The key 
to these theories is to depth-integrate general three-dimensi- 
onal equations in the shallow direction, resulting in a system 
of two-dimensional equations which still retains some infor- 
mation about variations in thickness. 

Let Oxyz be a coordinate system with the x-axis down- 
slope and the z-axis normal to a channel with mean slope 
Q. For simplicity, we further consider boundaries, flows, and 
external forcing to be (statistically) uniform in y. The contin- 
uum macro-scale fields are thus independent of y, while the 
DPM simulations remain three-dimensional and will be pe- 
riodic in y. The free-surface and base location are z = s{x,t) 
and z — b{x), respectively. The thickness of the flow is thus 
h{x,t) — s{x,t) — b{x), and the bulk density and velocity 
components are p and u = (m,v = 0,w)', respectively, as 
functions of x,y,z and t. 

The three-dimensional flow viewed as continuum is de- 
scribed by the mass and momentum balance equations ( fTST i 
and (fT9] l. At the top and bottom surface, kinetic boundary 
conditions are satisfied: D{z~ s) / Dt = and D{z~ b) / Dt = 
at their respective surfaces, and with material time deriva- 
tive 

D{-)/Dt^d{-)/dt + ud{-)/dx + wd{-)ldz 

(since we assumed v — 0). Furthermore, the top surface is 
traction-free, while the traction at the basal surface is essen- 
tially Coulomb-like. We decompose the traction t = t, + ?„ n 
in tangential and normal components, with normal traction 
tn — — n • (crn), where n is the outward normal at the fixed 
base. The Coulomb Ansatz implies that t, = — /x|t„|u/|u| 
with friction factor jj. > 0. Note that ji generally can be a 
function of the local depth and the flow velocity. Its deter- 
mination is essential to find a closed system of shallow-layer 
equations. 

We consider flows that are shallow, such that a typi- 
cal aspect ratio e between flow height and length, normal 
and alongslope velocity, or normal and alongslope varia- 
tions in basal topography, is small, of order ^(e). Further- 
more, the typical friction factor /i is small enough to satisfy 
H = ff(e^) with 7 G (0, 1). We foflow the derivation of the 
depth-averaged swallow layer equations for granular flow 
by 0GTNO3I without assuming that the flow is incompress- 
ible. Instead we start the asymptotic analysis from the di- 
mensionless form of the mass and momentum conservation 
equations (fTSl l and ( fT9] l. assuming only that the density is 
independent of depth at leading order. Density, velocity, and 
stress are depth averaged as follows 

O^fOdz. (33) 

In the end, one retains the normal stress ratio K — Oxx/^zz^ 
the shape factor a = u^/if', and the friction jx as unknowns. 



The goal is to investigate whether these unknowns can be 
expressed as either constants or functions of the remaining 
shallow flow variables, to leading order in 0'{e). The latter 
variables are the flow thickness h — h{x,t) and the depth- 
averaged velocity S = ii{x,t). At leading order, the momen- 
tum equation normal to the base yields that the downward 
normal stress is lithostatic, o;z(z) = pgcos9{s — z) + ^(e). 
Depth-averaging the remaining equations, while retaining 
only terms of order (?{epL), yields the dimensional depth- 
averaged shallow-layer equations, cf., ||VATK+07|IBT12L 

^ {hpu) + ^ (^hpau^ + ^gh^P cos 9^ = ghpS, (34b) 
with 

u db 

S = sinO - 11 — cosO - - — cos0. (34c) 

\u\ dx* 

To demarcate the dimensional time and spatial scales, we 
have used starred coordinates. These scales differ from the 
ones used before in the particle dynamics and the dimen- 
sionless ones used later in the DPM simulations. The shallow- 
layer equations ( |34] | consist of the continuity equation ( I34al l 
and the downslope momentum equation ( I34bb . The system 
arises also via a straightforward control volume analysis of 
a column of granular material viewed as continuum from 
base to the free surface, using Reynold-stress averaging and 
a leading order closure with depth averages. 

While the mean density p can be modelled as a system 
variable by considering the energy balance equation, we will 
assume that it can be expressed as a function of height and 
velocity p{h,ii). Thus, the closure to equations (|34] | is de- 
termined when we can find the functions p{h,u), K{h,u), 
a{h,u), and jU(/i,m). In Section l972l we will analyze if and 
when DPM simulations can determine these functions. 



4.2 Granular friction laws for a rough basal surface 

The friction coefficient, ji, was originally IIHM84II taken to 
be a simple Coulomb type jj. = tan 5, where 5 is a fixed fric- 
tion angle. Note that in steady state for a flat base with b = 
0, the shallow-layer momentum equation ( I34bl i then yields 
fi = tan0. Pure Coulomb friction implies that there is only 
one inclination, = 5, at which steady flow of constant 
height and flow velocity exists. That turns out to be unre- 
alistic. Three parameterizations for jj. have been proposed in 
the literature. 

Firstly, Forterre and Pouliquen IIFP03I found steady flow 
in laboratory investigations for a range of inclinations con- 
cering flow over rough basal surfaces. They measured the 
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thickness hstop of stationary material, left behind when a 
flowing layer was brought to rest, with the following fit 



hstopjO) _ tan(g2)~tan(e) 
Ad tan(0) -tan(5i) ' 



di<e <d2, 



(35) 



where 5i is the minimum angle required for flow, §2 the 
maximum angle at which steady uniform flow is possible, 
d the particle diameter, and A a characteristic dimensionless 
length scale over which the friction varies. Note that hstop 
diverges for 9 — 5i and is zero for 9 — 82- For h > hstop, 
steady flow exists in which the Froude number, the aspect 
ratio between flow speed and surface gravity-wave speed 
(F = u/y/g cos 9h), is a linear function of the height. 



i-stop 



(0) 



for for 5i < < ^2, 



(36) 



where /3 and 7 are constants independent of chute inclina- 
tion and particle size. Provided one uses the steady state as- 
sumption pL — tan 9 to hold (approximately) in the dynamic 
case as well, one can combine it with (l35l l and (l36l l to find 
an improved empirical friction law 



II = li*{h,F) =tan(5i) 



tan(52) — tan(5i) 
l5h/{Ad{F + Y)) + l' 



(37) 



It is a closure for fi in terms of the flow variables, which 
has been shown its practical value. Note that in the limit 
5i — >^ ^2 = 5 the Coulomb model is recovered. 

Secondly, in an earlier version ||Pou99l , another, expo- 
nential fitting was proposed for hstop, as follows 



h'stopW 



A'd 



= ln 



tan(0)-tan(5i') 
tan(52) -tan(5() 



for di<9< 82 (38) 



with the same limiting behaviour, and primes used to denote 
the difference in the fit. It yields the friction factor 



-tan5[)e^^'''(''+y)i . (39) 



/i = li'{h,F) = tan d[ + (tan 62 - 

Equation ( |35] | did, however, prove to be a better fit to exper- 
iments and is computationally easier to evaluate. 

Finally, an additional correction to ( l36l l was made in 
IIJen06J to include dependence on the inclination accounting 
for enduring particle contacts in very dissipative frictional 
flows. 



tan2(0) 



hs,op{9) tan2(5i 



(40) 



for which we can use any appropriate fit for hstop- It leads 
subsequently to a more complicated evaluation of the fric- 
tion law for IX. We omit further details. We will compare our 
DPM simulations against these rules and fits for the rough 
basal surface, and extend it to smoother surfaces in the up- 
coming sections. 




Fig. 2: DPM simulation for approximated height H — 17.5, 
inclination 9 = 24° and the diameter ratio of free and fixed 
particles, A = 1, at time t = 2000; gravity direction g as indi- 
cated. The domain is periodic in x- and y-directions. In the z- 
direction, fixed particles (black) form a rough base while the 
surface is unconstrained. Colours indicate speed, increasing 
from blue via green to orange. 



5 Simulation description 

In this section, DPM simulations are used to simulate monodis- 
persed granular flows. Parameters have been nondimension- 
alised such that the flow particle diameter d =1, mass m=l 
and the magnitude of gravity g = 1- The normal spring and 
damping constants are A^" = 2 • 10^ and 7" = 50; thus the 
contact duration is tc = 0.005 and the coefficient of resti- 
tution is e = 0.88. The tangential spring and damping con- 
stants are k' — {2/7)k" and 7' = 7", such that the frequency 
of normal and tangential contact oscillation and the normal 
and tangential dissipation are equal. The microscopic fric- 
tion coefficient was taken to be /i' = 1 /2. 

The interaction parameters are chosen as in IISEG^Oll to 
simulate glass particles of 0. 1 mm size; this corresponds to a 
dimensional time scale of \fdfg = 3 . 1 ms and dimensional 
velocity scale \fdg = 0.03 1 ms^ ' . The above parameters are 
identical to the simulations of Silbert et ai, except that dis- 
sipation in tangential direction, /, was added to damp rota- 
tional degrees of freedom in arresting flow. Adding of such 
tangential damping removes all vibrational energy for flows 
otherwise arrested. The differences in the simulations with 
7' — y" reported here are minor relative to the case with 
7' = 0. Silbert et al. also investigated the sensitivity of the 
results to the particle interaction parameters tc, e, the ra- 



10 



— — — fl = 60' 




— — — fl = 50° 

— — — fl = 40° 

— — — fl = 30° 



500 1000 1500 2000 



t 

Fig. 3: Shown is the ratio of kinetic to mean elastic energy 
over time for H — 20, basal roughness X — I, and varying 
chute angles Q. Flow stops for inclinations Q < 20.5°, re- 
mains steady for 21° < < 29° and accelerates for 9 > 30° 
(dashed lines). 

tio k" /k', and jU' ; they found that while the density of the 
bulk material is not sensitive to these interaction parame- 
ters, the flow velocity increased with decreasing friction ji^'. 
Nonetheless, the qualitative behaviour of the velocity pro- 
files did not change. 

The chute is periodic and of size 20 x 10 in the x- and 
y-directions and has a layer of fixed particles as a base. The 
bottom particles are monodispersed with (nondimensional) 
diameter A. Various basal roughnesses are investigated by 
taking A = 0, 1/2 to 2 in turn, with A = as flat base. This 
bottom particle layer is obtained by performing a simula- 
tion on a horizontal, smooth-bottom chute. It is filled with 
a randomly distributed set of particles and we simulate un- 
til a static layer about 12 particles thick is produced. Then 
the layers of particles at height z G [9.3, 11]A are selected, 
remaining particles deleted, and the selected ones moved 
downwards by 11 A. The layer is thick enough to ensure that 
no flowing particles can fall through the base. Their posi- 
tions are fixed. 

Defining a 'filling height' H, an integer amount of about 
200 H particles are inserted into the chute. To insert a par- 
ticle, a random location {x,y,z) G [0,20] x [0,10] x [0,/] is 
chosen, where I = H. If the particle at this position over- 
laps other particles, the insertion is rejected, and the inser- 
tion domain is enlarged by increasing / to 7 + 0.01 to en- 
sure that there is enough space for all particles. The initial 
packing fraction is about p/pp = 0.3. Thus the particles ini- 
tially compact to an approximated height H, giving the chute 
enough kinetic energy to initialize flow. Dimensionless time 
is integrated between t £ [0,2000] to allow the system to 
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Fig. 4: Overview of DPMs for A = 1, with markers denot- 
ing the flow state at f = 2000: arrested •, steady o, and ac- 
celerating * flows. Grey dash-dotted lines mark thickness li 
for fixed H from H = 10,20,30,40 upward. The demarca- 
tion line is fitted to hstop{d) in equation (l35T l (solid line) and 
^'stopi^) ill (dotted line). Error bars mark intervals es- 
tablishing the demarcation line. Red crosses denote the de- 
marcation between arrested and accelerating flow found in 
MSEG+Om . 

reach a steady state. A screen shot of the system in steady 
state is given in Fig.|2] 

To ensure that the size of the periodic box does not in- 
fluence the result, we compared density and velocity profiles 
of the flow at an angle = 24° and height H — 17.5 for do- 
main sizes Lv ^ 10, 20, 40, Ly = 10 and = 10, 20, 40, 
Ly = Lx/2, and found no significant changes. 

6 Arrested, steady, and accelerating flow 

From the experiments of Pouliquen IIPou99L granular flow 
over a rough base is known to exist for a range of heights and 
inclinations. DPMs by lISEG+Oll also showed that steady 
flows arose for a range of flow heights and (depth-averaged) 
velocities or Froude numbers. Their simulations did, how- 
ever, provide relatively few data points near the boundary 
of arrested and steady flow to allow a more adequate fit 
of the stopping height. In this section, we therefore per- 
form numerous DPMs at heights and angles near the sep- 
aratrix between the steady flow regime and the regime with 
static piles. To study the full range of steady flow regimes, 
simulations were performed for inclinations Q varying be- 
tween 20° and 60° and approximated or 'filling heights' 
H = 10, 20, 30 and 40. In Section[8l we will repeat (some 
of) these simulations for varying base roughness. 

Whether a steady flow regime has been reached, is quan- 
tified here by assessing the time whereafter the ratio of ki- 
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netic energy normalized by the mean elastic potential en- 
ergy becomes time independent. This is shown in Fig. [3] 
where we plot such an energy ratio for a rough base, con- 
stant height, and varying chute angle. The elastic potential 
energy is averaged over t G [1000,2000] to minimize fluctu- 
ations after start-up, but any interval larger than 100 appears 
sufficient. For chute angles at most 20.5° the kinetic energy 
vanishes after a short time, thus the flow arrests; for chute 
angles between 21° — 29°, a constant value is reached, indi- 
cating steady flow; and, for inclinations above 29° the en- 
ergy keeps increasing: thus flow steadily accelerates. If the 
energy ratio remained constant within t E [1800,2000], the 
flow was deemed steady, otherwise the flow was deemed to 
be either accelerating or stopping. 

Unlike fluids, the free surface of granular flows, and thus 
the flow height, are not well defined. In I SEG^Oll . the height 
of the flow was estimated by H, which is equivalent to as- 
suming a constant packing fraction of p/Pp = 7z/6. How- 
ever, the exact height h — s — b of the flow varies from the 
approximated height H due to compaction of the flow and 
H is typically an overestimate. In llVATK+07l . the surface 
of the flow was defined by the time-average of the maxi- 
mum vertical position of all flow particles. One could also 
define the free surface of the flow as the height where the 
density vanishes. The latter two methods, however, have the 
disadvantage that saltating particles can lead to slightly over- 
estimated flow heights. 

Instead, we will define the height via the downward nor- 
mal stress. For steady uniform flows the downward normal 
stress is lithostatic, i.e., balances the gravitational weight, 
such that 

(yzz{z) = p{z!)gcosedz'. (41) 

This is a direct consequence of the momentum balance equa- 
tions. Thus, C7--(z) has to decrease monotonically; the base 
and free surface are the heights at which G^ziz) reaches its 
maximum and minimum value, respectively. However, in or- 
der to avoid effects of coarse graining or single particles near 
the boundary, we cut off the stress Ozziz) on either boundary 
by defining threshold heights 

zi = minjz : cr,- < (1 — K)ma\a,A and (42a) 

^" zgR 

Z2 = max{z : Gzz > ifmaxcj^} (42b) 

with K=\%. We subsequently linearly extrapolate the stress 
profile in the interval (zi,Z2) to define the base b and sur- 
face height s as the height at which the linear extrapoation 
reaches the maximum and minimum values of a,,, respec- 
tively, 

b = zi--^—^iz2-zi), i = Z2 + Y3^fe-zi)- (42c) 



The variable most sensitive to these height choices is p. It 
shows well-defined functional behaviour for our definition 
of height, shown later. This is not the case if we define height 
by the density or the mefliod in llVATK+07l . The flireshold 
K = 1% was chosen because the results in Fig. fT2l were rel- 
atively insensitive to the choice of K at or above 1%. 

To determine the demarcation line hs{9;X) between ar- 
rested and steady flow with good accuracy, we performed a 
set of simulations with initial conditions determined by the 
following algorithm. Starting from an initial 'filling height' 
H = 4 and inclination 9 = 21°, the angle was increased 
in steps of 1° until eventually a flowing state was reached. 
From this initial flowing state, the height was increased by 
2 particle heights, if the flow arrested, or else the angle de- 
creased by 1/2°, assuming the curve is monotonically de- 
creasing. Flow was defined to be arrested when the energy 
ratio Ej^in/Egia fell below 10^^ within 500 time units, oth- 
erwise the flow was classified as flowing. In simulations in 
which such arrested flows were continued after t — 500, a 
further decrease of kinetic energy was observed, thus vali- 
dating the approach. This procedure yields intervals of the 
inclination angle for each height and, vice versa, height in- 
tervals for each angle, between which the demarcation Une 
lies. The values presented in |SEG+01 1 deviate at most 0.5° 
from our observations, perhaps due to the preparation of the 
chute bottom, or the slightly different dissipation used. A 
demarcating curve between steady and arrested flow was fit- 
ted to equations ( l35b and ( |38] | by minimizing the horizon- 
tal, respectively vertical, distance of the fit to these intervals, 
see Fig.m Fitting hstop{9) yields better results than h'^iapiO) 
for all roughnesses and only the fit ( |35] | will be used here- 
after Similar fits will be made in Section[8]for varying basal 
roughness. It leads us to a study of the depth profiles for 
steady state flow in the following section. 

7 Statistics for uniform steady flow 

To obtain detailed information about steady flows, we use 
the statistics defined in Sect. [3] Since the flows of interest are 
steady and uniform in x and y, density, velocity and stress 
wiU be averaged over x, y and f . The resulting depth pro- 
files will depend strongly on the coarse-graining width w, 
which needs to be carefully selected. Representative depth 
profiles for particular heights, inclinations and basal rough- 
nesses will also be analyzed. 

Depth profiles for steady uniform flow are averaged us- 
ing a coarse graining width w over x e (0,20], y G (0, 10] 
and t e [2000,2000 + r]. The profile of a variable x is thus 
defined as 

{X)liz)^7^^ / / Xwit,x,y,z)dxdydt, 

2007 J2000 Jo Jo 

(43) 
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Fig. 5: Depth-averaged norm of the momentum rate of 
change, r = \ d,{pu)\2dz, with d,{pu) determined by (l44l i 
for varying time averaging interval T. Steady flow at height 
H — 20 and incHnation = 24° was used. Temporal fluc- 
tuations decrease inversely proportional to the length of the 
time averaging interval. 



with Xw in turn the macroscopic field(s) of density, momen- 
tum and stress, as defined in (fTSl l. ( fT6b and (l32b . We average 
in time with time snapshots taken every ft /2 units. 

To determine an appropriate time averaging interval T, 
we calculate the rate of change in momentum from the den- 
sity, velocity and stress fields by 



djpu) 
dt 



= V • cr - pg - u • V(pu) 



(44) 



For steady flow, the temporal variations in mass and mo- 
mentum should approach zero when averaged over a long 
enough time interval T. This is shown in Fig. |5] where we 
plot the depth-averaged norm of the momentum rate of change 
for varying time averaging interval. For T > 100, the tempo- 
ral fluctuations decrease to less than 2% of the largest term, 
pg, in the momentum equation. In the remainder, we choose 
T — 100 as the averaging interval. 

The effect of varying coarse-graining width w is shown 
in Fig.|6l which shows the z-profile of particle volume frac- 
tion p /pp, where Pp is the particle density. For small w we 
observe strong oscillations of about 0.9 particle diameters 
width, particularly at the base. The microscopic oscillations 
are increasingly smoothed out and finally vanish as we ap- 
proach w — 0.5. For larger w, such as w > 1, the macroscopic 
gradients at the base and surface are smoothed out, an un- 
wanted effect of the coarse-graining. The same behaviour is 
observed in the stress and velocity fields. Smoothing over 
the microscopic effects makes it impossible to observe mi- 
croscopic layering in the density, which we still wish to 
identify in our averaged fields. Hence, we choose w = 0.25 
as the coarse-graining width, such that layering effects re- 
main visible along with macroscopic gradients. 

The microscopic oscillations at the base indicate a strong 
layering effect of particles near the boundary, despite the 
rough bottom surface. The macroscopic density throughout 
the flow is almost constant in the bulk and decreases slightly 
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Fig. 6: Particle volume fraction p /pp for H — 30, 9 — 24°, 
and A = 1 for varying coarse graining widths w. Micro- 
scopic layering effects are visible for w < 0.5. The density is 
constant in the bulk and decreases slightly in the basal layer. 
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Fig. 7: Normal and shear stresses for H = 30, 9 = 28°, and 
A = 1 . Shear a^z and downward normal stress a^z are bal- 
anced by gravitational forces, see equation iAli . The other 
normal stresses show anisotropic behaviour 



at the base. An approximately constant density profile is a 
feature of all steady flows and is an important assumption of 
depth-averaging. 

Non-zero stress components are plotted in Fig. |7] We 
observe that the stress components are nearly symmetric. 
Shear stresses Oyx and a,,., are zero since the flow veloc- 
ity in 37-direction vanishes. For steady flow, the downward 
normal stress Gzz{z) is lithostatic and satisfies equation (l4Tl l 
with a maximum error of 0.5%. Since the density is nearly 
constant, we obtain a linear stress profile, another assump- 
tion of depth-averaged theory. Applying the momentum bal- 
ance ST% to steady uniform flow further yields that the shear 
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Fig. 8: Flow velocity profile of thick flow for H = 30, 
=24°, A =0,2/3,5/6, 1,2 and fori/ = 30, =22°, 26°, 
A = 1/2. For a rough base with A > 5/6, we see a Bagnold 
velocity profile (dashed line), except near the surface. For a 
smooth bases with X < 2/3, the profile becomes more con- 
vex. For X — I /2, 9 < 24°, the flow velocity shows layering 
while still observing the Bagnold profile. For A = 0, a con- 
siderable slip velocity is observed. For A = 2, the basal shear 
is small due to flow particles trapped between basal particles 
so that the definition of the base b{x) is rather fuzzy. 



stress satisfies Gxz — jr p{z')g&^n9 dz' . Thus, the macro- 
scale friction jj. satisfies fi = (7^^/(7^-. = —gx/gz — tan 0. This 
relation is locally satisfied for all steady flow cases to an 
accuracy of |0-tan-i(ju)| < 0.4°. The remaining normal 
stress components, Gxx and (7„,, are not constrained by this 
mass balance. We thus see in Fig. Q significant anisotropy 
in the amplitude of the normal stresses, in particular in ff„., 
showing that the confining stress is largest in the flow direc- 
tion, except for very small inclinations. It is always weak- 
est in the lateral or y-direction with fluctuations at the base 
that are in phase with the fluctuations of the density. Gen- 
erally, the anisotropy increases with higher inclinations and 
smoother bases, as will be analyzed in future work. 



8 Transition from rough to smootli base 

Next, we study the effect of smoother bases on the range of 
steady flows by decreasing the diameter A of the base parti- 
cles, with the limiting case of a flat bottom wall for A = 0. 
Such an extensive numerical study of the effects of chang- 
ing bottom roughness appears to be novel. To that effect, the 
DPM simulations from Section |6] were extended such that 
results for basal roughnesses A =0, 1/2, 2/3, 5/6, 1, and 2 
can be compared. Before we show the /is,op-curves for these 
simulations, we investigate the extent to which changes in 
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Fig. 9: Demarcation line hs{9;X) for varying basal rough- 
ness. Markers denote the midpoint of the intervals around 
which the curve was fitted. Steady flow is observed at 
smaller inclinations for smoother bases. While the smaller 
angle 5^ i varies only slightly, the larger angle 82 f^, de- 
creases rapidly with the smoothness. For A = 0, the demar- 
cation line is vertical at = 12.5°. 



basal roughness lead to more complex density and velocity 
profiles. 

We summarize the density profiles seen without explic- 
itly showing the results. For decreasing basal roughness A, 
we observe that the microscopic oscillations and the dip in 
density at the base increase, while the bulk density remains 
constant. For A = 1/2 and A = 2/3 and small inclinations, 
we see steady flow that is strongly layered throughout the 
flow. In contrast, for A = 2 there is a low flow density in 
the basal region, since some of the free particles are small 
enough to sink a little into the base, forming a mixed layer 
of fixed and free particles. 

Velocity profiles for H — 30 and = 24° are shown for 
varying basal roughness in Fig. [8] For A = 1, we observe 
the Bagnold profile (Bag54] for thick collisional flows, dif- 
fering only at the surface. For very thin flows at H ~ 10 or 
inclinations near the arresting flow regime, the profile dif- 
fers strongly from the Bagnold profile and becomes linear. 
For smoother bases, the flow velocity increases, and the pro- 
file becomes more concave. Weak to stronger slip velocities 
are observed for A < 2/3. For A = 0, thicker flows have con- 
stant velocity throughout the depth, almost corresponding to 
plug flow. 

A family of demarcation curves hstop ( 9 ; A ) between steady 
and arrested flow is shown in Fig.|9] The curve fits are based 
on 



h„op{9;X) =Axd 



tan(52,A) -tan(0) 
tan(0)-tan(5i,;i)^ 



(45) 
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in which the dependencies on A are expHcitly denoted. The 
fitting parameters 5i x, 82^1, Ax appearing in (05]) are found 
in Table [T] Again, a fit based on the original equation ( |35] | 
(or (l45T l) rather than Pouliquen's early fit dJSl ) yields the best 
results. 

For a flat bottom, such that A = 0, steady flow initiates 
and resides at or very tightly around an inclination = 12.5° 
for all heights, see Fig. |9] It is in agreement with the angle 
found in the laboratory experiments of IIGTDD03I . Hence, 
for a smooth base the flow is steady only at a single inclina- 
tion, arrests for lower inclinations and accelerates for larger 
inclinations. Such behaviour is in line with laboratory obser- 
vations and DPM simulations, e.g. [VATKtOlJ- For 1/2 < 
A < 2, we observe Pouliquen-style behaviour in Fig.|9] The 
angle 5i i of flow initiation is nearly constant with respect to 
X. In contrast, the range of angles at which both steady and 
arrested flow is possible, 82,1 ~ ^i.A- is maximal for A = 1 
and decreases for smoother chutes with A < 1, as follows 
from Table [1] This has been reported in IIGTDD03II for lab- 
oratory experiments, who also observed a slight decrease of 
the interval 82^1 — 5i x for A > Ac « 2. However, their A^ 
was measured for basal particles fixed at the same height 
and depended on the compactness of the base. We observe a 
slight decrease of 82^1 for A = 2; however, the fitting curves 
in Fig.|9]do mildly overlap for A > 1. 

We recall that di x and 82,1 are fitting parameters for the 
hstop-cnrve ( |45] ) which does not necessarily imply, though it 
is expected, that the flow accelerates for angles greater than 
82 X ■ Surprisingly, while steady flow is observed exclusively 
for € (5i,i,52.i) when A = 1, the range of angles associ- 
ated with steady flow for smoother chutes (i.e., when A < 1) 
extends to greater inclinations with 6 > 52,x ■ For these latter 
cases, 5„cc i > 82^x is defined as the smallest angle at which 
accelerating flow is observed; the DPM simulations show 
fliat 

/25°±1° ifA=0, 

Oacc 1 = ■{ (46) 

[29° ±1° otherwise. 



A 




S2.1 


Ax 




Ta 


err 





12.25 


12.25 




1.500 


-4.065 


0.886 


1/2 


17.913 


21.357 


11.959 


0.300 


-0.236 


0.229 


2/3 


17.217 


24.302 


11.088 


0.215 


-0.143 


0.104 


5/6 


17.425 


26.828 


7.465 


0.203 


0.039 


0.113 


1 


17.708 


32.780 


3.355 


0.196 


0.040 


0.121 


2 


17.518 


29.712 


5.290 


0.189 


0.080 


0.137 



Table 1: Table of fitting parameters 5^ x, 82 x^ for the 
curve hstop{0',X) and Px, 7a for the flow rule (l47b . including 
the variance of the flow rule, err{F — F^ata), for all steady 
flows. 
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Fig. 10: Top: Overview of DPM simulations for A = 1/2, 
with markers denoting the flow state at f = 2000: arrested •, 
layered x, oscillating o, steady o, and accelerating * flows. 
Demarcation line hstop{d'A f^) is fitted according to ( |35] |. 
Bottom left panel: Profile of particle volume fraction of lay- 
ered flow at // = 20, = 22°. Bottom right panel: Ratio 
of kinetic over mean elastic energy for oscillating flow at 
H = 30, 6 = 24°. An example of oscillating flow is seen 
Fig.S 

Note that for the A = case, between angles 12.5° and 25°, 
the flow is steady and layered, because the friction factor is 
nonzero. 

The above is illustrated in Fig. [TO] for A = 1/2, where 
one observes the two different steady state regimes. At higher 
angles, ^3 1/2 < < 5acc, 1/2^ a disordered regime similar to 
that for a rough base is observed. At smaller angles, 5i 1/2 < 
G < ^3 1/2, the flowing system self-organizes into a state of 
layered flow consisting of ordering in the jc-y-plane for the 
bulk (bottom left panel of Fig. [TOl i. except for a small inter- 
mediate region, 6 ~ ^3 1/2, where a transitional flow regime 
can be found. It is characterized by large oscillations in the 
ratio of bulk averaged kinetic to elastic energy due to a spon- 
taneous ordering and disordering, or stop-and-go flow, of the 
system as a function of time (lower right panel). The same 
flow regimes have been observed in IISGPL02I . where the 
smoother bottoms were achieved by arranging the base par- 
ticles in a grid-like fashion. In contrast, we always use of a 
fully disorded base and vary the roughness by changing the 
basal particle size. 
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9 Closure relations for the depth-averaged model 

The goal of this section is to close the shallow-layer equa- 
tions ( |34] | by a determination of the basal friction /i, the 
mean density p, the stress ratio K, and the velocity profile 
a, using our DPMs. A demarcation will be made of the flow 
regimes in which such a determination is possible. 



9.1 Friction pL of shallow-layer model 

For the rough base several friction laws have been proposed, 
as detailed in Section l4!2l In the following, we will compare 
these friction laws for the rough base, ^ = 1, as well as for 
varying basal ratios X. 

To obtain a function for the basal friction pL, we used the 
approach of PouUquen, who found that for a rough base the 
Froude number is a linear function of h/hstop{0)- A first ap- 
proach was to fit the Froude number to a linear function of 
h/hstop{0\X) across the range of non-accelerating DPMs. 
While this does work for A > 5/6, a (linear or other) fit 
does not work for well for A < 2/3 because for the smoother 
bases layered and oscllating flows occur for 5^ i < Ostop {h,X) 
< dj x - This is illustrated for A = 1/2 in Fig.[TOl Instead, 
the Froude number is fitted with hstop{0'A) such that 



h,topiOA) 



-7a, for 5iA<0<5a. 



23° ±1° ifA^O, 

25.5° ±0.5° if A 1/2 and// = 10, 

53jt 24.5°±0.5° ifA = l/2and//>10, 

24.5°±0.5° if A =2/3 and// = 10, 

, Ostop {h\X) otherwise. 



error for the fit to WT\ are found in Table [T] with a standard 
error defined by 



err(Mti) = (£x2/(A^-l)) 



(47) 



where 5^ x is the largest angle at which oscillating flow is 
observed. In Fig. [10] these angles 5^i^^.x and 5t, x are shown 
for the case A = 1/2. Overall, the simulations reveal that 



(48) 



1/2 



(49) 



We remark that a fit to equation ( l36b is marginally better than 
Jenkins' adaption (|40] |. but the differences are too small to 
discriminate accurately. 

The situation for layered and oscillating flows is more 
complicated. We illustrate that for the case A = 1 /2. Two fits 
are shown in Fig. [TT] one for the layered case (dotted line 
concerning the crosses), and one for the steady case (solid 
line concerning the circles). The oscillating flows seem to 
defy a sensible fit because the flow swings irregularly be- 
tween the layered and disordered states. That oscillating be- 
havior was also shown in Fig. [TOl (bottom right panel). 

For steady flow, the shallow-layer equations (l34b yield 
pL — tan ((J). Indeed, within the range of steady or arrested 
flows DPMs confirm directly that the friction at the base lies 
within \\L — tan^^((J)| < 0.4°. In summary, for the steady 
flow regimes observed in our DPM simulations, the friction 
< coefficient of the depth-averaged equations (|34|) is parame- 
terized to be 



p.{h,F;X) =tan(5i.i) 



tan(52,i) — tan(5] 



i,ij 



l5xh/{Aid{F + n)) + V 
for d,^x<d<SaccA- (50) 



where the parameters 82.1, Ai are independent of the 
base; and, Px and yx are depending explicitly on A . All val- 
ues are found in Table [T] with Sacc.A. and 5^ x given in ( |46] | 
and (l48T l. Despite its determination for steady flows, such 
a closure for jj. is assumed and often observed to be a rea- 
sonable 'leading order' approach for unsteady shallow-layer 
flows. Furthermore, for smoother bases, closure laws for lay- 
ered and oscillating flows have eluded us. It seems that the 
homogenization and steadiness assumptions of depth-averaged 
shallow-layer flow break down in these cases. 



The results of such fits to the Pouliquen law for 8^ x < d < 
dacT.X ^6 shown in Fig. (TT] with corresponding fitting pa- 
rameters provided in Table [T] Shown is the Froude num- 9.2 Functions p,K, a of shallow-layer model 
ber F = S/^/gcos Oh against the ratio of flow and stopping 



heights h/hstopiO; I). For the disordered steady flow regime, 
concerning angles 53^^ < < 5„£.(Jl, the data are seen to 
fit better with the stopping angle hstop{0\X = 1), the one 
for basal surface A = 1, rather than with the actual stop- 
ping height hstup{0',X). This is a key observation. It shows 
that the Froude number F increases as the roughness A de- 
creases, due to the lower dissipation at the base. The weaker 
Froude number dependence for A = 2 seen in the right panel 
of Fig.[TT]is in line with the zero shear observed at the base 
in Fig. [8] The full set of fitting parameters and the standard 



DPM simulations of steady uniform flows are considered for 
disordered steady flow with 5^ < 9 < 5^^^., to determine 
closures for p, K and a as functions of continuum fields 
u and h. The layered and oscillating flow regimes are thus 
excluded momentarily. 

All steady disordered flows show a constant density pro- 
file in the bulk of the flow, cf. Fig. |6] while the density de- 
creases near the base and the surface. The lower density 
region at the base spans about two particle diameters for 
A > 0, and less than 9d for A = 0, while the surface region 
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V/!.„op(e,A = i) h/h,,,p{e,x = i) h/Ku,p{e,x = i) 

Fig. 11: Froude number/^ = u/^fgh over height scaled by the stopping height for /L = 1 (left), A = 1/2 (centre), and for 
all basal roughnesses (right). Data with symbols 'x' denote steady layered, 'o' steady, and 'o' oscillating flows. Data with 
symbols '.' correspond to arrested or steady flows near hstop- The data is fit using h^f^pi^Q ^ \ — 1) (solid lines). 



spans always less than Ad. Thus, a mean bulk density can 
roughly be defined as 

h-bd Jb+2d 

In Fig. [12] the bulk volume fraction and the mean volume 
fraction are shown for roughness A = 1 and varying height 
and inclination. The bulk volume fraction decreases with in- 
clination 0, but is independent of flow height and rough- 
ness, whereas the mean volume fraction depends also on 
flow height and roughness. We fit the mean bulk density of 
all steady disordered flows with A > to an arbitrary func- 
tion 

P/" IPp = Co - exp ( - c2 ) /c 1 , (52a) 
with fitting parameters 

CO = 0.5985, ci =4.8°, andc2== 40.85°. (52b) 

Standard deviations of the mean bulk volume fraction and 
mean volume fraction for all cases with A > are 

eiT(p/''-pc) =0.002, and err(p/'' -p) = 0.018. (52c) 

Secondly, the normal stress ratio K — Oxx/^zz is deter- 
mined. It describes the anisotropy of the stress tensor and is 
expected to be unity under isotropic conditions. The range 
of K for steady disordered flow is generally small, ranging 
from 0.98 to 1 .07, except for A = 0, where it can be as low 
as 0.68. The stress anisotropy generally increases with incli- 
nation. For X >Q,K fits to a function linear in 0, 

Kf" = \ + {e-dx)/do (53) 

with do = 132° and di — 21.50°. The model results give a 
small standard error of eniK-K^'') = 0.013. Given that 
the dependence on inclination is small, we may as well take 
KKi 1. 

Finally, we develop a fit for the shape factor a(A) = 
iP'/u^. The fit is based on a phenomenological model of the 
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Fig. 12: Mean volume fraction p/Pp for roughness A = 1, 
and varying approximated heights H and inclinations 9. The 
mean volume fraction in the bulk, p / Pp, denoted by *, col- 
lapses onto a function of the inclination (solid line), while 
it shows a small dependence on the flow height, due to the 
density decrease near base and surface. 



observed velocity profiles, as shown in Fig. [8] For rough 
bases A > 5/6, a Bagnold velocity profile. 



mb(z) = ^« (l - {ih-z)/hf^^ , (54a) 



is observed in the bulk of the flow; a linear profile in the 
surface layer, which is about 5d thick; and, a convex profile 
with no slip in the base layer, whose thickness increases 
as the height approaches the stopping height. No kinks occur 
at the intersection of the layers. Thus, we model the velocity 
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Fig. 13: Left figure: fitting parameter as a function of hstop{0; 1 ) //i for X ~ I and varying height h and inclination 9. The 
solid line shows the linear fit used to obtain a from equations ( l54l i. Right: Shape factor a for A = 1 and varying height h and 
inclination 9. Markers denote the simulation data, while dotted lines denote fits using i54\ with corresponding coefficients 
from Table 12] Fitted values and simulation data are connected by a solid line. 
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Fig. 14: Depth profile of normalized strain, {h/u)d^u cor- 
responding to velocity profiles shown in Fig. [8] For rough 
bases, the strain is modelled by a Bagnold profile, except 
near the base and surface. For smoother bases, X < 2/3, the 
profile becomes linear. For A = 0, a large slip velocity is ob- 
served, and the strain becomes inversely proportional to the 
depth. 

by 

^(z;feA)-< ^(z), b^<z<m?ixis~5,b 

[ ^max{s — 5,bx), otherwise, 
M(0;/7A)=0forA>5/6. (54b) 

A fit to the strain d^u in Fig.[T4]shows it is well approximated 
by this model. The parameter b;^ decreases with increasing 



distance from the stopping height, and a simple fit reads 

bx^d^h,rop{9,l)/h, (54c) 

where hstop{d,l) was chosen since hstop{9,X) does not pro- 
vide values for all inclinations for which steady flow is ob- 
served. Subsequently, the shape factor a(A) = m^/ can be 
computed numerically and compared to the measured values 
in Fig. [13] The coefficients b^ are given in Table |2] 

For A < 2/3, the dependence of the shape factor on height 
and inclination diminishes and can be approximated with a 
constant value a (A). The Bagnold profile disappears and 
the flow becomes more convex and plug-like, as shown in 
Fig. [8] Each velocity profile will be analyzed in turn next. 

For A = 0, the slip is so large that we can assume plug 
flow to hold. There is almost no shp for X —2/3 and a slip 
of approximately u{Q) / max^u{z) = 1/6 for A = 1/2. We 
neglect the variations at the surface and the bulk and model. 
Thus, we model the velocity profiles as 

(5/3^{l-z/h)\ A =2/3, 

^ = <^ 0.16 + 0.84(5/3 -(l-z//;)2), A = 1/2, (55) 
[l, A = 0, 

The corresponding coefficients a (A) are found in Table |2] 
and provide a good fit to the data. 

In summary, the functions p, a and K depend on the in- 
clination 9 and the height h. The inclination 9 in turn can 
be written as a function of the friction coefficient ji such 
' that 9 = tan^^(/x(/!,M)). This allows us to describe the pa- 
rameters of the shallow-layer model in terms of the height 
h, roughness A, and friction jj, {h, u) and thus provides a clo- 
sure proper for the system. The different behavior for the 
varying A's remains an open issue, since we only provided 
emperical fits above. 
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A 






err 





1.000 




0.012 


1/2 


1.142 




0.01919 


2/3 


1.201 




0.02741 


5/6 




5.105 


0.02558 


1 




6.678 


0.01763 


2 




13.75 


0.04479 



Table 2: Fitting for the shape factor a = a (A) for A < 2/3 
and a = iP'/u^, u — u[z;bx) for A > 5/6, and the standard 
error. Closure relations are fitted to all data sets of steady 
unordered flow, 5^ < 9 < 5^^.^. 

10 Conclusion 

10.1 Summary 

In this article, an extensive series of DPM simulations was 
used to determine closure relations for the popular shallow- 
layer model of granular flows on inclined chutes. The latter 
model is a depth-averaged continuum model with a macro- 
scale variable thickness h = h{x,t) and with a mean velocity 
u = u{x,t) as variables. For simplicity, we assumed unifor- 
mity in the lateral y-direction. The flow consisted of monodis- 
persed particles of diameter d and the base of monodispersed 
particles of diameter At/. The bottom roughness, or diameter 
ratio, A was systematically varied. Simulations revealed the 
existence of a range of chute inclinations 6 for which hori- 
zontal and temporal variations are small enough to produce 
an approximately steady and uniform flow. Particle flows 
with variations in height h and inclination 6 were numeri- 
cally investigated for varying basal roughness A . 

We observed the following phenomenology: at small in- 
clinations, the flow quickly forms a static pile, while at large 
inchnations, the flow continues to accelerate. Between these 
two regimes there was a range of inclinations in which steady 
flows occurred (cf., Fig.|4|i. The curve hstop{0',^), a function 
of height versus inclination, forms the separatrix between 
arrested and steady flow, as a function of basal roughness 
(Figs. |4] and |9]l- For smaller basal roughness, steady states 
arise at smaller inclinations and heights, and the range of an- 
gles shrinks for which steady flow is possible. Other types 
of steady flow were observed at small inclinations for small 
base particles, showing a strong layering in depth as well as 
disordered and oscillatory flows {cf., Fig.llOb. 

Depth profiles for density, velocity and stress were con- 
structed using coarse-grained macroscopic fields. The coarse- 
graining width was carefully chosen to preserve some mi- 
croscopic structure as well as macroscopic gradients (cf., 
Fig.|6]l. The assumptions of depth-averaged theory were con- 
firmed in the simulations for a certain range of steady, uni- 
form flows: the density was constant at depth, and the down- 
ward normal stress as well as shear stress were lithostatic. A 
often-used key assumption, often used, is that statistically 



steady DPM simulations, or laboratory experiments, are rel- 
evant to find the closure relations even for time-dependent 
continuum shallow-layer models. 

Consequently, four closure relations could in principle 
be determined: for basal friction jj., stress ratio K, mean den- 
sity p, and for the shape of the velocity profile a. Firstly, 
basal friction = tan 9 was shown to be a function of height 
and flow velocity (Table[Tll. Pouliquen's approach was found 
to be valid, with the Froude number as a linear function 
of h/hstop{Q',^ — !)■ This fitting approach was extended 
to smoother cases with A < 1, where the Froude number 
was fitted to h/hsiop{9;X = 1) instead of the h/hsiop{9;X) 
for the actual A or basal roughness. The stopping curve as- 
sociated with the diameter A = 1 of the flowing particles 
is more relevant than the stopping with the actual A. One 
possible explanation is that there is a boundary layer of in- 
termittently slow flow particles that originated in the bulk, 
and that shields the smoother base from the bulk flow. Clo- 
sure relations for the mean density p, stress anisotropy K 
and shape factor a were also established as follows. For 
rough bases with A > 5/6, the determined closures were 
vaUd for Ostop{h;X) = 5^,1 < 9 < 5„C6-,A. with 9stop{h;X) 
the inverse of the hstop{9; X)-cnrve between arrested and 
dynamic flow. For smaller roughnesses with A < 2/3 and 
9stop{h;X) < 9 < layered and oscillating flows arose 
for which we are (as yet) unable to capture closures. For 
these smoother bases with A < 5/6, closures were obtained 
for the range <S < 5„cc,A- 

10.2 Open questions 

What does the granular shallow-layer model enable us to 
do, and what can we not do with it? In the range of steady 
flows, this continuum model can be used to predict steady 
and time-dependent flows. Strictly speaking, this is only al- 
lowed for steady flow in the established inclination range 
^3. A < 9 < Siic^.i, but it can be expected to remain vaUd 
for the slowly-varying dynamic cases as well. It is often the 
case, however, that even rapidly-changing flows can be cap- 
tured by models that should only be valid for the slowly- 
varying cases. Consequently, a systematic study of the va- 
lidity of the shallow-layer model is required. By respec- 
tively extending the "hydraulic" analysis for fluidized granu- 
lar matter and water in Vreman et al. ||VATK+07|| and Akers 
and Bokhove IIAB08I . granular flows within constrictions 
become a nice, analytic ally- treatable targets. Such flows in 
constrictions reach a steady state and appear (partially) ac- 
cessible by direct DPM simulations. 

What do these results enable us to do? Whether the steady 
DPM-based closures are valid across granular "hydraulic" 
jumps in such steady and constricting flows is of interest. 
Whether the steady DPM-based closures hold for (slow) tran- 
sient routes towards such steady states is of interest, too. 
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What closures should be used outside the formal range of 
applicability for the smoother bases, so for the layered and 
oscillating flows for Ostopih;X) < 9 < 5^ i and the acceler- 
ating flows for 9 > appears a tantalising, and as yet, 
open question. 

What are we not able do? Although, we did observe 
layered and oscillating flows in our DPM simulations, it is 
doubtful as to whether the homogenization assumption that 
led to the shallow-layer model is sufficient. Nonetheless, the 
lithostatic balance relation is shown to hold for the DPM 
simulations, as expected from standard asymptotic analysis 
using the aspect ratio of normal to planar velocity and length 
scales. 



10.3 Outlook 

Alternatively, a multi-scale modelling approach might be 
adopted such as the heterogeneuous, multiscale methodol- 
ogy rtWEL+07l . among others, in which closure relations 
for discretisations {e.g., [?]) depth-averaged shallow-layer 
models are coupled to DPM simulations in selected regions 
in space and time. Thus computational costs would be di- 
minished while accurate closure relations are gathered inter- 
mittently in time and space. 

For future work, we advocate the extension of our DPM 
simulations with investigation of the three-dimensional clo- 
sure relations. We surmise that reduced lithostatic models 
for shallow granular flows could be more consistently de- 
rived from three-dimensional continuum models with stress 
closure determined from DPM simulations in combination 
with laboratory measurements. These new models would be 
reduced and therefore computationally still manageable for 
large-scale debris flows; for example, when the degrees of 
freedom in the vertical remain limited, but are extended be- 
yond only one degree of freedom. Such reduced modelling 
is akin to hydrostatic modelling in water-wave and coastal 
hydrodynamics. 
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A System of Hamiltonian equations in tlie 
dissipation-free limit 

The purpose of this appendix is to show that the particle system in 
the dissipation- and yield-free limit, i.e., 7" = 7' = and ^1^ — )■ 0°, 
is a Hamiltonian system. It subsequently facilitates the derivation of 
so-called conservative or symplectic discretization schemes, in time. 
These have been shown and are believed to provide better long-term 
statistics than classical time discretization schemes. Furthermore, anal- 
ysis of the Hamiltonian limit allows one to clearly demarcate the trans- 
fer of energy between its kinetic, elastic, and internal components. 

We show that the normal and tangential forces are elastic, that is 
the system does not dissipate energy; instead kinetic energy is con- 
verted into potential energy in the springs and vice versa. If the tan- 
gential spring is not fully unloaded when two particles loose contact, 
the potential energy stored in the tangential spring is converted into 
internal potential energy in each particle (vibrations). 

We use the notation given in ij2] In the dissipation- and yield-free 
limit, the contact force between particles /, j is given by 



(56) 



where 5," and SJy are given by equations {T) and {9)- The equations of 
motion for translational and angular momentum of particle i are given 

by. 



d , d 



dt' 



(57a) 
(57b) 
(57c) 



Mi 



in three dimensiones, where r, is the position, ct, is the angle, p, the 
momentum and 0, the angular momentum of particle (. 

To define the Hamiltonian system, we pair these generalized posi- 
tion and momentum vectors as follows 



Q(0 = {r,.(0,«,-(0}f=., P(0 = {P,'(0,<^,'(0}f=i- 



(58) 



Then the kinetic energy can be calculated using only the general- 
ized momenta P as follows 



(59) 



The potential energy is a combination of the potential of gravity, the 
potential of the normal and tangential springs and the internal poten- 
tial energy in the particles, created from the remaining potential of the 
tangential spring at the time tfj that a particle pair j} looses contact. 



(60) 



with r,y = ry(f) — r,{f). The potential can be expressed in terms of the 
position and the tangential springs at all times, which itself is a function 
of the previous positions of the particle pair, 

V(Q) = VgraAQ)+VeU,{Q)+Vi„,{Q). (61a) 

where the gravitational, elastic and internally stored potential energy is 
defined by 

N 

V,v(Q) = XI g (61b) 



1=1 

;=l j=i+i \ ^ ^ 
N N /,! 

Vi„m = Y. I I {-T\^Wn 



(61c) 
(61d) 



We will now show that the total energy H = T + V satisfies the Hamil 
tonian equations, 

dH _ dpi dH _ dcpi 
(9r, dt ' dcxj dt ' 



dH dti ^ dH doLi 
dpi dt ' dcpi dt 



(62a) 
(62b) 



To derive j62at . we calculate 



dH__d_ 
dti dti 



(63) 



term by term. We can show that 

d^~2 " max(0,(J; + rfy)/2 -'■;.,■) 



= rmax(0, {di + di)/2-rij) 



d((di+dj)/2-rij) 



dn 



= -k" max (0, {di + dj)/2 - nj) 

''ij 

= -l<"5Finij, 



(64) 

where we used the fact that max (0, (rf, + dj)/2 — r/j) is continuous in 
time. 

Further, we take the derivative d,..^5'jj (a denotes the vector coor- 
dinate) by the chain rule and use ^21' to show that 



d5':j _ dt dS\^ 
dria{t) dru,{t) dt 

1 



Via^a + (-Viaea ■ n//)n// + (5 ■ ■ • ( - V/„ e,, ) ) 



51 '*'./ 



(65) 



where denotes the a-th basis vector of the coordinate system. Here, 
i65) becomes 



d k< 
dr:,, 2 



dS< 



5''=k'S'r- 



(66) 

where the cancellation of terms arises because the tangential spring is 
orthogonal to the normal vector. After substituting i66) and J64t into 
l l63t we obtain 

dH d ^ ^ k",.„,2 k\., ,2, 

_ = _(_,„,r,g+^_|5/)|^ + -|5;,|^). 



= -mrg-'^k'^SI'jnii-k'S'ij = 

Mi 



dPi 
dt ■ 



Next, we calculate 



dUi,, dau, 



Mi 



Mi 



(67) 



(68) 



We take the derivative using the chain rule and equations (9) 

and O 



d8< 



dUiait) dai„{t) dt 



dt dS]j 



dt d{-COie„xbij) 



dUiait) 
-e„ X by. 



dt 



(69) 
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Substituting (69) into j68J we obtain 

rlH 

^ = -^^%.(.„xb,,) = i:(^'5;,xb,,.)„, 

where we used the identity 

c- (e„ X d) = -(c X d)„ Vc,d e 

Thus, using 468b and that b,y and n,j are parallel, we obtain 

d<i>i 



B Orthogonality of the tangential spring 

(70) To show that the tangential spring is orthogonal to r,y, note that the 
tangential spring is initially of zero length and therefore orthogonal to 



(71) A^'ij-'iii) 



dt 



dt 



dH 
dcti 



Subsequently, we derive 462b t since 

dH _ d Ipip _ Pi _ dtj 
dpi dpj 2m, ntj dt ' 

and 

dH d (pi dcxi 

depi ^ dcpi 21, ^ li ^ dt 



dt 



(72) 



(73) 



(74) 



■Tij + S'ij-yij 



(77) 



Finally, we show that the total energy is conserved. Since mass 
m,, radius r, and spring constants k", k' are constant, H has no direct 
dependence on / and thus 



Thus, we can integrate equation ill) to obtain a continuously orthogo- 
nal tangential spring with 6'. ■ r,j = 0. 



C Algorithms for time integration and the calculation of 
the tangential force 

The algorithm for the time integration and the calculation of the tan- 
gential force is shown in Algorithms[T]and|2] 



dH_ 

IT 



0. 



(75) 



Using and J621 and M5\ yields 

d , , dH ^ dH dti dH da, dH dp, dH d<pi 
dt dt A'j dvi dt doij dt dpi dt d<pi dt 

dH dpi dti dcf>j dcti dti dp, dec, d<^, 
dt dt dt dt dt dt dt dt dt 
= 0, (76) 

Fig. ^] shows the energy balance of two particles colliding non- 
coUinear in the dissipation- and yield-free case. One can see the jump 
in energy at the end of contact, where potential tangential spring energy 
is converted into internal energy. 
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2 I- 
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Algorithm 1: Time integration 

Data: Initial positions and translational and angular velocities 

rj', v-', ajj*, masses »!,■, inertias /,■, time step At 
ti ^ ri, V,- ^ v?, a;,- ^ a;f V/ 
(f,,q,) i — forces-and-torques({rj,Vj,ajj}^^[) V/ 
for;= l,2,...,A'do 



V/ i — 
(jji ^ 
ti i — r,' + Atyj Mi 

foreach particle pair (;, j) in contact do 



if contact is 



new then 



■ 5;.+4?a;, 



(f„q,) 

V/ i — 



forces-and-torques ( {ry , y j , t^; j ] ) V/ 

2 m: 



y^ + ^^Mi 



u^i + ffVi 




Algorithm 2: Calculation of the tangential force, in- 
cluding sliding 







0.2 0.4 0.6 O.i 



1.2 1.4 1.6 



r..^-i^si.-Yy'.. 

if(|l^,.|>Mc|17yl)then 



Fig. 15: Energy balance of two particles colliding non- 
collinear in the dissipation- and yield-free case. The kinetic 
and elastic potential energy are in balance, until the parti- 
cles loose contact and potential spring energy is converted 
into internal energy. 



